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The purpose of the research project "Multi-Disciplinary Optimization of 
Aeroservoelastic Systems" was to advance the methodology of multi-disciplinary 
optimization and to develop efficient analytical and computational tools for 
simultaneous structural and control system optimal design with aeroservoelastic 
constraints. These tasks have been accomplished successfully. The new 
analytical methods, computation codes and case studies developed during the 
course of this work extend the knowledge of this complex subject and supply the 
aircraft developer with practical guidelines and efficient computational tools 
for achieving integrated design goals in an optimal way. 

The first part of the work, in which the Minimum-State aeroservoelastic 
modeling method has been extended to deal with large problems with varying 
parameters, has been reported in the interim report of March 1st, 1990. The 
main results of this part are summarized in Appendix A, which has been accepted 
for publication in the AIAA Journal, and utilized in the MIST computer code. 

The second part of the research project included further developments of 
analytical and computational tools and a synthesis of an end-to-end optimization 
procedure from data-base construction to the numerical optimal design process. 
This part is summarized in Appendix B, which will be presented at the 3rd 
Airforce/NASA symposium on multi-disciplinary analysis and optimization in 
September 24-26, 1990, San Francisco, CA, and utilized in the FASER computer 
code. 
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Extensions to the Minimum-State Aeroelastic Modeling Method 

Mordechay Karpel* 
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Introduction 

> 

In order to account for unsteady aerodynam.es in first-order, time-invariant 
state-space formulation of aeroelastic equations of motion, the aerodynamic forces have to 
be described as a rational function in the Laplace domain. The Minimum-State (MS) 
aerodynamic approximation method 1 - 2 was designed to minimize the number of aerodynamic 
states in the resulting aeroelastic model. References 2 and 3 applied the MS method to 
subsonic aeroservoelastic problems with one flutter mechanism and demonstrated a reduction 
of about 75% of the number of aerodynamic states relative to other methods with the same 
level of accuracy. The effectiveness of the MS method was increased by the introduction of 
a physical weighting technique 2 which weights each aerodynamic input data term according 
to its relative importance. Reference 4 used the MS formulation for additional reduction 
of the model size by dynamic residualization of high frequency structural states. The MS 
and the physical weighting procedures are extended in this paper to expand their 
efficiency and generality and to improve the dynamic residualization. Even though the 
formulation and numerical examples deal with structural-mode-related aerodynamics only, 
the extensions are applicable to control surface and gust related aerodynamics as well. 


The MS Approximation Procedure 

The MS method approximates the Laplace domain generalized aerodynamic foice 
coefficient matrix by: 

[Q s (p)J = [A 0 ] + :Aj]p + [A^P 2 + [D](p[I] - [R]) _1 [E]p (1) 
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where p is the nondimensionalized Laplace variable p=sb/V, where b is a reference length 
and V is the true airspeed. The resulting time-domain state-space aeroelastic equations of 
motion are presented in Ref. 2. The number of aerodynamic states (m) is equal to the 


order to [R] - 

The input data are unsteady aerodynamic complex matrices [FCk^l+itGCk^)], 

calculated at several p=ik l points where each k=(DbfW is a tabulated reduced frequency. 
The approximation problem is to find the combination of the real valued [A Q ], [AJ, [A^, 
[R], [D] and [E] of Eq. (1) that best fit the tabulated data. The mxm aerodynamic lag 
matrix [R] is diagonal with distinct negative values to be chosen by the analyst. The 
applications of Refs. 2 and 4 indicated that the results are not very sensitive to the lag 
values when they are spread over the range of tabulated ^ values. Three approximation 
constraints are applied to each term of [Q $ ] in order to reduce the problem size by 
explicitly determining [A Q ], [Aj] and [A 2 ]. The formulation of Ref. 2 is extended here to 
allow more flexibility in constraint selection without increasing the problem size. The 


three constraints are: a) data match at k-0, which yields: 


V. - F ij (0) 

1J 


( 2 ) 


b) real-part data-match at a non-zero k =k r or a zero coefficient constraint, which yield: 

‘■'’ 1 (3) 


A 2.. 

ij 




F--(O) - FjAfil/kJ + lDj| T fk|[I]+[R] 2 


IE.) or A, = 0 

J lj 


where {D-} T is the ith row of [D] and (Ed is the jth column of [E]; and c) imaginary-part 
1 J 

data-match at a non-zero k x =k g , or a zero coefficient constraint which yield: 


A. =G..(kJ/k +|D.| T fki[I] + [R] 21 

hi 11 g 


-1 


[R](E-1 


or 


A l.. = ° 
d 


(4) 


'ij'“g" g ' l B 

The approximation formula (1) and the constraint Eqs. (2), (3) and (4) yield an over 
determined set of approximate equations which are solved for [D] and [E] by an iterative, 
weighted least-square procedure which starts with an initial guess of [D]. The equations 
and the solution procedure are those of Ref. 2, modified to allow k f of Eq. (3) and k g of 
Eq. (4) to have different values for different aerodynamic terms, and to allow the data 


2 



match constraints to be replaced by zero coefficient constraints. 


Approximation Constraints for Subsequent Dynamic Residualization 
The MS aeroelastic model is used in Ref. 4 for a further reduction of the model stze 
via dynamic residualization which eliminates the stales associated with a subset of high 
frequency v.bratton modes, but retains most of their effects on the retained states. The 
coefficient matrices of Eq. (1) are partitioned into the retained (r) and eliminated (e) 

partitions: 


i t A i) 


rA. 

A. 1 

for i = 0,1,2 ; [D] - 

[D 1 
r 

i 

rr 

1 

re 

D 

A. 

i 

L er 

A. 


e J 

x ee-l 




; [El = (E E e ) 


(5) 


Unlike the static residualizat.on, which neglects all the e related partitions except for 
the Aq terms, the dynamic residualization neglects only the A^. A^, and A^ 

terms The retamed effects of Aj , A, , D £ and E e .improve the accuracy of the 

er re 

residualized model without increasing its size. The attempt made in Ref. 4 to improve the 

dynamic residualization even further by constraining the neglected terms to be zero ,n the 

preceding MS procedure (in lieu of data match constraints) did not yield better results. 

The reason was that with A, =0 the approximation errors are increased significantly. 

ee 

This had a negative effect on the quality of the enure approximated aerodynamics 
(including that of the retained modes) because the MS procedure minimized a single total 
error parameter. The modification suggested here is to apply the least-square solutions 
for [E ] and [D ] with the data associated with the retained modes only. The [D e ) and 
(E ) matrices are solved with the entire data. As a result, the approximated [Q S J is 

not affected by the inclusion of the eliminated modes in the approximation procedure. 

Physical Weighting 

A physical weighting method which weights each term of the tabulated aerodynamic data 
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according to a "measure of importance” is presented in Ref. 2. The measure-of-importance 
matrix associated with [Q t .(ikJ] is: 

v) V 

•v 

[W] = f-[MJk^ + i[B s ]k t + [K $ ] + q n tQ s (ik i )] 

where [M ], [B,] and [K $ ] are the generalized mass, damping and stiffness matrices and q n 

is a nominal dynamic pressure. As shown in Refs. 2 and 3, the variations of the 

measure-of-importance terms of Eq. (6) with k may have very sharp peaks. In addition, the 

peak values of many terms may be several orders-of-magnitude smaller than other peaks. 

The resulting extreme variations of weights may cause unrealistic approximation curves. 

To ensure good results at k values which fall between the tabulated ones, and to 

facilitate the application of the resulting aeroelastic model to a variety of flow 

conditions, structural modifications and control parameters, it may be desired to widen 

the weight peaks and to scale up the extremely low weights. The peak widening is 

performed in n wd cycles where, in each cycle, W^) is changed to be maxlW^), 

W..(k ), W..(k .)) of the previous cycle. The weights to be applied in the MS procedure 

\J ^ J 

are then calculated by: 





max 


W O 
1 cut ) 

maxtW.Y W.. J 
i , j J 



where W.. - max| Q (ik^) W..1 
q i c ij J J 

and where W is defined by the analyst. In this way, the maximum weighted absolute 
c u t 

value of each aerodynamic term falls between W cu[ and 1. 


Numerical Examples 

The numerical examples deal with the mathematical model of the Active Flexible Wing 

(AFW) wind-tunnel model (described in Ref. 4) with symmetric boundary conditions at Mach 

0.9. The Doublet Lattice tabulated oscillatory aerodynamic matrices were generated at 12 

k values between 0.0 and 2.0 using the STABCAR 5 computer code which was also employed to 
i 
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calculate the baseline p-plane roots using the p-k method with ten vibration modes. The 
resulting root-locus plots are shown in Fig. 1 which indicates two flutter mechanisms. 
The flutter dynamic pressure and frequency of the first mechanism (second branch) are 
1.447 psi and C0p58.94 rad/sec (k^ ut =0.21). The flutter results of the second 
mechanism (seventh branch) are qp4.247 psi; 0)^194.43 rad/sec (kfl ut =°-69)- 

The physical weightings were performed with q=1.2 psi. Two types of physical 
weightings are compared below. The first type (symbolized by P-0) is with the original 
measures-of-importance of Eq. (6), namely with no peak widening (n wd =0) and with no 
upscaling (W =0 in Eq. (7)). The second type (symbolized by P-2) is with the two 
peak-widening cycles (n wd =2) and with =0.01. The maximal weighted magnitudes of the 
P-0 aerodynamic data terms are given in Table 1. The most important modes are 2, 3, 6, 7 
and 10 which have the highest diagonal values. The off-diagonal values associated with 
these modes are also higher than those of most other terms. It can be noticed that about 
50% of the terms in Table 1 are smaller than 0.01. These terms are scaledup to 0.01 in 
the P-2 case. Comparison of Table 1 with the aeroelastic behaviour of Fig. 1 indicates 
that the weighting is reasonable over the entire q lange of interest. 


i= 

1 


r 

J. 

) 

3 

\ 

4 

5 


6 

7 

8 


9 

10 

i—1 

o 

098 

0 

.126 

0. 

086 

0.012 

0. 

000 

0.018 

0.007 

0 . 

001 

0.006 

0.005 

J 1 

2 

o 

137 

1 

000 

0 . 

707 

0.038 

0. 

007 

0.129 

0.104 

0 . 

008 

0.079 

0.045 

3 

0 

085 

0' 

.722 

0 . 

622 

0.052 

0. 

003 

0.091 

0.039 

0 . 

003 

0.029 

0.015 

4 

o 

010 

0 

.059 

0 . 

044 

0.022 

0. 

001 

0.006 

0.003 

0 . 

000 

0.002 

0.000 

5 

o 

001 

0 

.003 

0. 

.004 

0.002 

0. 

043 

0.002 

0.001 

0 . 

000 

0.001 

0.001 

6 

o 

027 

0 

.132 

0. 

. 115 

0.013 

0. 

002 

0.377 

0.026 

0 . 

003 

0.012 

0.060 

7 

0. 

013 

0 

.07 9 

0. 

. 053 

0.006 

0. 

001 

0.042 

0.274 

0 . 

004 

0.028 

0.066 

8 

o 

001 

0 

.008 

0 

. 005 

0.000 

0. 

000 

0.003 

0.003 

0 . 

018 

0.002 

0.001 

9 

o 

Oil 

0 

.070 

0 

. 041 

0.006 

0. 

.001 

0.023 

0.031 

0 . 

003 

0.135 

0.090 

10 

0. 

007 

0 

.046 

0 

.021 

0.007 

0. 

.001 

0.087 

0.087 

0 . 

000 

0.098 

0.291 


Table 1: Maximal weighted magnitudes of the aerodynamic data terms 


The flutter characteristics of the resulting state-space models have been found by a 
linear root-locus analysis with variable q. The quality of the approximations is evaluated 
by comparing the state-space results with STABCAR result. An overall measure in each case 
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is the RMS value of the percentage errors in the four flutter parameters ^q^ ,(0f ,q^ 311(1 

CO ). Comparisons between RMS flutter errors in non-weighted (N) cases and physically 
f 2 

weighted (P-0 and P-2) cases are shown in Fig. 2. It can be observed that the P-2 cases 
generally yield the best results and that they are more consistent than the P-0 cases. 

Calculations at Mach 1.15 (not shown) exhibit similar results. 

The root locus of the P-2 case with 6 aerodynamic lags, 

diag[R]=(-0.2,-0.45,-0.8,-1 2,-1.7,-2.0) is compared in Fig. 1 to that of the reference 
STABCAR solution. It can be observed that the agreement is good over the entire ranges of 
frequency, damping and dynamic pressure. This indicates that the physical weights 

calculated at q are adequate over the entire range. 

All the MS cases above were constrained to match the data at k=0.0 and at the highest 
tabulated reduced frequency, namely k^k g =2.0. Data match constraints at ^ values close 
to either one of the two k flut values caused a slight improvement in the respective 
flutter mechanism, but a slight increase of the overall RMS error measure by about 1%. An 
error increase of about 2% was obtained when [A 2 )=[0] replaced the k f constraints. Much 
more significant errors resulted from the replacement of the k constraints by ApO even 
when applied to the aerodynamic terms associated with the highest trequency mode only. 
These errors were reduced considerably with the application of the new procedure for 

subsequent dynamic residualization. 

To demonstrate the application of the modified MS approximations in subsequent 
flutter analysis with dynamic residualization, the P-2 model with 6 aerodynamic states has 
been extended to include the first 20 vibration modes (instead of 10). The MS 
approximations were performed with the special residualization constraints assigned to the 
last 10 modes. The reference case is flutter analysis with all the 20 vibration modes. 
Reduced-size flutter analyses were performed by eliminating a subset of high-frequency 
modes by either mode truncation, static residualization or dynamic residualization. 
Variations of flutter dynamic pressure percentage errors vs. number of eliminated modes 
are shown in Fig. 3 for the second flutter mechanism. Similar trends, but with smaller 
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errors, were obtained for the first flutter mechanism. These results demonstrate that MS 
aerodynamic approximations facilitate additional high accuracy model size reduction via 

dynamic residualization. 
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Figures 


Figure 1: Root locus of the 10 mode subsonic AFW case. 

K.ure 2: RMS flutter errors resulting from MS aerodynamic approximations. 

Figure 3: Flutter dynamic pressure errors of the second flutter mechanism vs. # 
of eliminated high-frequency modes. 
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Multidisciplinary Optimization of Aeroservoelastic Systems 


Mordechay Karpel* 

Technion - Israel Institute of Technology 

Haifa, Israel 


Abstract 

Efficient analytical and computational tools for simultaneous optimal 
design of the structural and control components of aeroservoelastic systems ar 
presented. The optimization objective is to achieve aircraft performance 
requirements and sufficient flutter and control stability margins with a minima 
weight penalty and without violating the design constraints. Analyti 
sensitivity derivatives facilitate an efficient optimization process whic 
Allows a relatively large number of design variables Standard finite-element 
and unsteady aerodynamic routines are used to construct a mo a 
Minimum-state aerodynamic approximations and dynamic residualization methods are 
used to construct a high accuracy, low order aeroservoelastic model. 
Sensitivity derivatives of flutter dynamic pressure, control stability margins 
and control effectiveness with respect to structural and control design 
variables are presented. The performance requirements are utilized by equality 
constraints which affect the sensitivity derivatives A gradient-based 

optimization algorithm is used to minimize an overall cost function, 
realistic numerical example of a composite wing with our con ro s 1 
demonstrate the modeling technique, the optimization process, an eir accuracy 

and efficiency. 


Introduction 


The design of the structural and control systems of a flight vehicle with a 
given aerodynamic configuration starts with separate analyses which are aimed 
satisfying basic stress and performance requirements. These preliminary designs 
are then optimized to satisfy structural and control, stability margin 
requirements with a minimal performance cost and without violating the design 
constraints. Modern high-performance, control-augmented aircraft may have a 
strong coupling between the structural and control systems through aeroelas ic 
effects. This calls for a multidisciplinary optimization process in which 
structural and control design variables are modified simultaneously. 


The purpose of this work is to present a practical and efficient 
optimization scheme in which the various aeroservoelastic aspects are analyzed 
and synthesized with a common model. The applicability of^ this approac o 
realistic design cases has been demonstrated by Livne et al. who presented an 
integrated synthesis scheme that can treat a rich variety of behaviour 
constraints and performance measures such as stress, displacements, contro 
surface travel and hinge moments, aeroservoelastic poles, gust response, drag 
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and aircraft maneuver parameters. A thorough literature survey is also given in 
Ref. 1. The structural model of Ref. 1 is based on the equivalent plate 
approach of Giles with which a complex structure can be represented by a 
relatively low number of degrees of freedom. This allows the inclusion of the 
entire structural model in the aeroservoelastic model, which is impractical with 
common finite-element models like those of NASTRAN, 

The optimization suggested in this paper can start with any structural and 
aerodynamic linear models, which has an advantage in practical applications. It 
is assumed that the structure can be represented in the optimization by a 
limited set of low-frequency vibration modes of the base-line model. The 
general scheme of a major optimization cycle is given in Fig. 1. An aeroelastic 
modal data base is first constructed for a relatively large number of modes. 

I. Data-Base 


*• 

t 



Fig. 1. The general scheme of a major optimization cycle. 
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The computational efficiency of this part is not as important as that of the 
later parts because it is performed once for the entire optimization process or 
at least for the major part of it. The aeroservoelastic modeling part of the 
scheme is aimed at resulting in a model of minimal size, but still accurate 
enough for the optimization process. Several size-reduction techniques are 
combined to achieve this goal and to provide the analyst with physical insig t 
and numerical measures for an a priori evaluation of the size-reduction effects. 
The third part of the optimization scheme is the actual model update process o 
minimize a cost function within the design constraints. High computational 
efficiency of this part facilitates an interactive design process where the 
design can inspect the results and change parameters in order to perform 
trade-off studies and to obtain realizable design. 

More details, analyticals development and numerical applications are given 
in the following sections. For the sake of clarity we start with the 
description of the test-case model that will be used to demonstrate the various 
parts of the optimization process. The formulation and.applications are limited 
in scope to allow elaboration on the new features of this work. 


The Test Case Model 

The test-case deals with the Active-Flexible-Wing (AFW) composite wind- 
tunnel model tested at NASA Langley Research Center. A top view of the NASTRAN 
finite-element model is given in Fig. 2. A description of the aerodynamic model 


y [in] 



x [in] 


Fig. 2. Top view of the AFW structural model. 


3 




1 

i 


appears in Ref. 3. The mathematical model assumes antisymmetric boundary 
conditions at Mach 0.9. The lowest 12 vibration frequencies and the associated 
mode descriptions appear in Table 1. 


i- 

Mode # 

Frequency (Hz) 

Description 

1 

0. 

Rigid body roll 

2 

7. 023 

1st fuselage bending 

3 

7. 856 

1st wing bending 

4 

13.069 

Missile pitch 

5 

16.161 

2nd fuselage bending 

6 

27.408 

2nd wing bending 

7 

38.271 

1st wing torsion 

8 

39.639 

Missile yaw 

9 

41.137 

Fuselage torsion 

10 

49.922 

3rd fuselage bending 

11 

51.640 

3rd wing bending 

12 

57.403 

2nd wing bending 

1 ---- 


Table 1: Natural Frequencies and Modes of the Base-Line Structure 

The subject for optimization is the wing-box structure between y=18.45 in 
and y=49.30 in. The wing-box composite upper and lower skins are mirror images 
of each other. The wing-box is divided into 11 optimization zones as shown in 
Fig. 3. The base line weight of the optimized portion of the structure is 3.06 
lb. The total number of high-strength graphite-epoxy plies in e o ach of o the skins 
varies between four at zone 2 (one in each of the 0 , +45 , -45 and 90 

o o o 

orientations) and twenty at zone 9 (four in 0 and two in each of the 28 , “62 , 
45°, -45°, 73°, “17° and 90° orientations). 

The model has 4 control surfaces per wing (see Fig. 2) driven by 
third-order actuators. A zero-order control system reads the output of a roll 
rate-gyro located near the centerline and commands the actuators through 
different gains. A performance analysis established that a control law which 
supplies a rolling moment of L=-3000 lb-in per unit aircraft roll velocity 
(minus pilot roll command) is required for adequate roll performance at the 
design dynamic pressure q^=1.5 psi. This yields the equality constraint 

T.G.7) .L = L/q, (1) 

l l 5 d 

i 

where L_ is the rigid aerodynamic rolling moment per unit dynamic pressure due 
o. 

1 

to unit deflection of the ith control surface and and T}^ are the associated 
control gain and aeroelastic effectiveness. In addition, the closed-loop 


4 





Fig. 3. Wing-box optimization zones. 


aeroservoelastic system is required to exhibit a sufficient flutter margin 

q f > 1-44 q d (2) 

All the non-zero gains; are required to exhibit sufficient gain margins 


GM.I >6 db uj 

l 1 _ 

where GM.=20 log(G?/G.) where G*/G. is a positive factor by which G has to be 
i ~ i i li 

multiplied to achieve instability at (with other gains remain unchanged). 

GVG >1 yields the positive gain margin and G*/G <1 yields the negative one. 
i i li 

Phase margin requirements and their sensitivity derivatives are given in Ref. 4. 

They are not discussed here because they are not critical in the test-case 

problem of this paper. 


The Aeroelastic Data-Base 

The aeroelastic data-base (see top part of Fig. 1) has been constructed 
using the NASTRAN code. The modes consist of the 25 lowest frequency vibration 
modes of the base-line structure, [\p )» and 4 control surface deflection modes, 
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f 
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[ip ] f of unit rotation of the respective control surface and zero deflective 
c 

elsewhere. The modes serve as generalized coordinates in the Laplace domain 
equation of motion 



2 

s + 


B 


s + 




Q (s) 
s 




2 

s 


+ 


r 




Q (s) 
c 



(4) 


where [M ] [B ] and [K ] are the generalized structural mass, damping and 
s 1 s s 

stiffness matrices respectively, [M ] is the control coupling mass matrix, {£} 

o 

is the' vector of generalized structural displacements, {5} is the vector of 

control surface commanded deflections, [Q (s)] and [Q (s)] are the generalized 

s c 

unsteady aerodynamic force coefficient matrices and q is the dynamic pressure. 

For the base-line structure [M ], [K 3 and [B ] are diagonal where the first two 

s s s 

were calculated from the finite element model and [B ] was constructed with 

structural damping of 0.01. [M ] can be calculated by 

c 



(5) 


where [M] is the mass matrix in discrete coordinates. A convenient way to 

calculate [0 ] and [M ] is by disconnecting the control surfaces from the 

c c 

actuators in the finite-element model and adding degrees-of-freedom representing 
{5}. These degree-of-freedom are mounted in the normal modes analysis by a 

fictitious mass matrix [M^] of large magnitudes. The test-case analysis was 

i 6 

performed with the 4x4 [M f ] matrix equal [ I] 10 while the maximum value of the 

other terms in the mass matrix is less than 1. The resulting modes include [$ ) 

as rigid-body modes while the other modes, frequencies and generalized masses 
remain practically identical to those of the original structure. Each control 
surface mode is normalized to a unit deflection of the respective <5^. 

Orthogonality implies 


</> 1 
m 

I 

" M 

0 


‘♦o' 

I 

Wh 

E 

_i 


0 

■» 

M i- . 


. I . 


[ 0 ] 


( 6 ) 


where the matrices are partitioned according to the original coordinates and the 
added ones. Equations (5) and (6) yield 



which is of course much more computationally efficient than Eq. (5). 


The unsteady aerodynamic code is employed to calculate [Q ] and [Q ] of Eq. 

iD L 

(4) (using [0 ] and [0 ] ) for various values of the nondimensionalized Laplace 
_ m c 

variable s=ik where k is the reduced frequency o»b/V where w is the vibration 
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frequency, b is a reference semi-chord and V is the flow velocity. The 
aerodynamic data-base of the test-case include 13 generalized aerodynamic 
matrices for the k^values of 0., 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 1.0, 1.4, 

1.9, 2.5 and 3.1. 


The remaining part of the aeroelastic data-base is the data needed for the 
structural optimization. The unit value of a structural design variable p g in 

our case represents one composite fiber ply added to the upper and lower 

surfaces of one of the zones of Fig. 3. Three design variables, representing 

the 0°, 45° and -45° orientations are assigned for each zone. The stiffness and 

mass matrices, [K. ] and [M. ] associated with p are calculated by the 

i l s i 

finite-element code and pre and post multiplied by the appropriate partitions of 

[\b ] to yield 
m 


d[K ] 

s 




d[H ] 
s 




( 8 ) 


Expressions for the derivative of [M ] with respect to p are given in Ref. 4. 

c i 

In our case 3[M ]/5p = 0 because the structure of the control surfaces is not 

c s. 

l 

modified throughout the optimization. The geometrical changes due to adding 
plies to the wing surface (stacking effects) are assumed to be negligible. 

Aeroservoelastic Modeling 

The aeroservoelastic modeling process (see part II of Fig. 1) is aimed at 

obtaining an efficient constant coefficient, state-space model to be used in the 

simultaneous structural and control optimization. To do so, the tabulated 

[Q (ik )] and [Q (ik )] matrices have to be approximated by rational functions, 
sc ci 

A review and extensions to the most common rational approximation methods and 
the associated model formulations are given in Ref. 5. Among those, the 
Minimum-State (MS) method 6 ' 7 has been shown in several applications to realistic 
problems 7 ' 8,9 to yield the most efficient aeroservoelastic models per desired 
accuracy. Best MS results are obtained when the physical weighting algorithm of 
Refs. 7 and 10, which weights the tabulated aerodynamic term according to their 
aeroservoelastic importance, is applied. The [Q (ik )] terms are weighted 

according to their effect on the determinant of the system matrix on the left 
hand side of Eq. (4), calculated at the tabulated s=ik b/V points, and the 

[Q (ik )] terms are weighted according to their effect on the return signal in a 

C l 

Nyquist analysis performed with all the gains having the same non-zero values. 
The aerodynamic matrices associated with the 25 vibration modes and the 4 
control modes were approximated with 6 MS approximation roots, which yield only 
6 aerodynamic augmenting states. 

The formulation of the resulting aeroservoelastic model is developed in 
Refs. 4 and 7. The closed-loop equation of motion for stability analysis reads 


7 



{x} = [A]{x> 


(9) 


where [A] is the closed-loop system matrix and {x} is the state vector combined 
of the structural states {£} and {£}, the aerodynamic augmenting states {x & } and 

the control system states {x >. In our case, the full-size model has n=68 

c 

states (50 structural, 6 aerodynamic and 12 control states). The eigenvalues of 
the constant coefficient, real-valued matrix [A] are used to analyze the system 
stability. Root locus analysis with variable q yields the flutter dynamic 
pressure, q , at which one of the root branches crosses' to the right hand side 

of the Laplace domain. The gain margins of Eq. (3) are found by a similar 
analysis with variable G. which yields the instability gain G*. More details 

and numerical examples for the calculation of flutter, gain and phase margins 

are given in Ref. 4. The accuracy of the MS approximation in the test-case has 
been examined by comparing flutter results with those of NASTRAN (using the PK 

method). The differences in q and w f were less than 1%. 

♦ 

' Model Size-Reduction 


The computational efficiency of the optimization process is approximately 
proportional to n . It is therefore desired to reduce n as much as accuracy 
allows. The vibration modes are divided for this purpose into three groups, 
those which may be truncated, those which may be eliminated via residualization, 
and those which remain in the model as independent states. The physical 
weighting algorithm described above and in Ref. 7 has been found to be a helpful 
guide in the mode selection-process. Two measures of aeroservoelastic importance 
are assigned to each mode. The measures are 


Q 



max (w. . 

j, ijt 


Q (ik ) 

S, . C 

1J 


} 


and 


( 10 ) 


Q* 

2 . 

l 


max 


r 


j. h 


w. . 

ijt 


Q 


.. ( Vl} 
1.1 J 


where W. . 

1JL 

tabulated 


is the physical weight 

[Q (ik ) ] or (Q (ik )]. 
sc c c 


assigned to the (i,j) 
The modal measures 


term of either the 
of aeroservoelastic 


importance of 17 of the 25 data-base modes are given in Table 2 in decreasing 

order of overall importance. The remaining 8 modes, for which both Q* and Q* 

i i 

are less than 0.01, are truncated in the optimization process. The 6 modes in 
the 2nd group of Table 2 are eliminated via residualization and the 11 modes in 
the first group of Table 2 are retained as independent variables. 
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i-: ::-—- i 

Mode 

t 

# 


Q* 

z 

Mode 

# 

q ; 

Q* 

2 

1 


0. 482 

1.000 

5 


0.022 

0. 035 

4 


1.000 

0. 125 

8 


0. 031 

0. 009 

3 


0. 611 

0. 404 

16 


0. 033 

0. 003 

6 


0. 512 

0. 035 

18 


0.024 

0. 007 

12 


0. 264 

0. 034 

24 

■j 

0. 023 

0. 005 

2 


0. 149 

0. 114 

25 


0.024 

0. 003 

11 


0. 161 

0. 035 

22 


0.019 

0. 004 

7 


0. 175 

0. 023 

13 


0. 017 

0. 004 

9 


0. 135 

0. 026 





_i 

t 

The 

other 

modes have 

Q* values of less 

than 

0. 01 



Table 2: Modal Measures of Aeroservoelastic Importance 


The truncation is performed by crossing out the associated rows on columns 
in [A] of Eq. (7). The truncated [A] satisfies 

( 11 ) 

the retained and residualized states 
respectively, and {U} is the complex column eigenvector related to s=io)^. The 

residualization in this work is performed by the dynamic method of Ref. 3. The 
residualized system matrix [A] satisfies 


A 11 A 12 
A 21 A 22 



f IF 1 

1 

1 

\ r = 1( V 



IF j 

J i 2 J 

^ 2 J 


where subscripts 1 and 2 relate to 



( 12 ) 


It is assumed that the differences between and {U^} of Eq. (11) and those of 

Eq. ,(12) are negligible. Consequently, the flutter eigenvalue and eigenvector 
of [A] are used to calculate {U^l from the bottom partition of Eq. (11). 


The effects of truncation and residualization on q^. and its derivative with 

respect to p are shown in Fig. 4. The solid lines give the percentage errors 

S 2 

introduced by truncation vs. the number of truncated modes. The base-line 

values of q =1.736 psi and dq /dp =0.188 psi have been obtained with the full 
* * S 2 

25 mode model where the gain values associated with the four control surfaces 
are 0.0, -0,1, “0.1 and 0.0664 respectively. The modes truncated at each point 
of Fig. 4 are those with the lowest measures of aeroservoelastic importance 
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Fig. 4. Size-reduction effects on flutter dynamic pressure. 

defined in Eq. (10). It can be deduced from the error rates of change that the 
usage of these measures for mode selection is adequate. The errors due to the 
truncation of 8 modes and residualization of 6 modes, also given in Fig. 4, 
demonstrate the accuracy of the dynamic residualization relative to that of 
truncation. 


Control Effectiveness 

The expressions for the aeroelastic control effectiveness parameters of Eq. 
(1) and their sensitivity derivatives follow the modal approach of Ref. 11. The 


real-valued [Q (0)], [Q (0)] and [K ] of Eq 

s c s 

and elastic mode partitions: 


Q (0) 



s’ 

; |q (0) 

s 

■ — 


is, 

%2 ■ 

L c J 


and 


(4) are partitioned into the rigid 



(13) 
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- - 

K 


" 0 

0 

s 


0 

mm 

K 

s J 


When there is a single rigid body mode in roll, the effectiveness of rolling 
moment due to unit deflection of the i-th control surface is 



where {Q } 0 
c 2 


is the j~th column of 


[Q c ] 2 


and 


(14) 


lB] = [K s l-[4K s ] + q d [C! s ] 22 , ^ 




a[K ] 

_ s 

= Z P s. 

IS. 1 


(15) 


where 3[K ]/3p is calculated by Eq. (6) using elastic modes only. The 
s s. 

1 


differentiation of Eq. (14) with respect to a structural design variable p 


s. 
1 


yields: 


dr). 


dp 


s. 

l 


Q 




" ' 

% 

12 

B 

mm m 

- - 


n aim ) 

-1 s 


dp 


i. 
j 


s. 
1 


B 


1 Q 


c[2 . 
J 


(16) 


It should be noted that the calculation of all the control effectiveness 
parameters in Eq. (14) and all their derivatives in Eq. (16), for a given set of 
design variables, involves only one matrix inversion of order n where n is the 

number of elastic modes after truncation (16 in our case). 


The variation of r). and dr)./dp , divided by their base-line values, with 

4 4 s 0 

the number of low frequency modes taken into account are shown in Fig. 5. The 
base-line values, calculated by all the data-base modes are v =0.216 and 


/3p =0.0079. 

a s 


The stars indicate the modes which will be actually truncated 


in the optimization. The convergence rates demonstrate the adequacy of the 
modal approach in considering control effectiveness in structural optimization. 
It should be noticed, however, that the convergence of the sensitivity 
derivatives is slower than that of tj itself. 


Sensitivity Derivatives and Constraints 

The sensitivity derivatives of flutter dynamic pressure and gain and phase 
margins with respect to structural and control design variables are presented in 
Ref. 4. They are based on the neutral stability eigenvalue, the associated 
column and row eigenvectors and derivatives of [A] of Eq. (7), namely the system 
matrix before residualization. When the model is residualized, the full-size 
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Fig. 5. Control effectiveness variations with number of structural modes taken 
into account. 


column eigenvector is calculated from the reduced-order one as discussed after 
Eq. (12). The full-size row eigenvectors are calculated in a similar way but 
with Eqs. (11) and (12) transposed. 

Equality constraints are defined as a linear dependency between the design 
variables of the form 


ZC.p. = C (17) 

11 o 

One of the design variables, p^ which has a non zero C^, is defined as a 

dependent variable. The differentiation of Eq. (17) with respect to an 

independent design variable, p., yields the constrained derivatives 





ac. 

C. + 2 -5-1 
l J 5 




(18) 


In our test-case, the rolling moment requirement of Eq. (1) is utilized by 
defining as a dependent gain, which yields the constrained sensitivity 
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derivatives with respect to other gains 


6 


dG . 
J 


a 


dG , 
J 


”j L 6 


j 3 


r 3G. 

4.L 4 

4 5 4 


( 19 ) 


and the constrained sensitivity derivatives with respect to the structural 
design variables 



dp 


s . 
J 







L a. G i 

1 



( 20 ) 


The Optimization Process 

Three cases of the base-line structure with different control gains are 
given in Table 3. All the cases satisfy the rolling requirement of Eq. (1) with 
G 1 =0.0. Case 1 is with G 2 =G 4 and G 3 =0 > which yields q f <q d - Since dq f /dG^ is 

relatively large and r\ is relatively low, can be used to suppress flutter 

while the other gains compensate for the loss of rolling moment. This is 
utilized in cases 2 and 3 in which q f >q d but the flutter and gain margin 

requirements of Eqs. (2) and (3) are not met. Even though the main participants 
in the flutter mechanisms of cases 2 and 3 are the same modes (3 and 4 of Table 
1), the flutter frequency, gain margins and the derivatives are substantially 
different, which indicates a high sensitivity of the stability characteristics 
to parametric changes. 


To obtain a balanced structural design, all the design variables^ associated 

with the -45° ply direction are constrained to be equal to the +45 ones. To 

avoid violation of stress requirements and to obtain a smooth design, all the 

structural variables are limited to 0<p <p where p =1.0 in zones 1 and 2 of 

s. s s 

l 


Fig. 3, 
control 


2.0 in zones 3 to 5, 3.0 
gains are constrained to 


in zones 6 to 8 and 4.0 in zones 9 to 11. The 


G =0 
1 


and the other G. <G where G is equal to 


either 0.09 or 0.1. The remaining independent design variables are the twenty 

o 


two 0 and 45 structural variables and two control gains (G, 


and G ). 


A preliminary analysis indicated that the critical stability parameters are 

q and the positive GM . values for j=2 to 4. Consequently, the cost function to 
f J 

be minimized is 


J 


A^AW-AVM V^fcrV 
e + e 



e 


(GM -GM.) 
o J 


j=2 


( 21 ) 


where q,. and GM are 
f o 

o 


the requested flutter dynamic pressure and gain margins, 
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f -- -J 

parameter 

Case 1 

Case 2 

Case 3 

G^ 

-0.0871 

-0.0900 

-0. 1 

2 




G^ 

0. 0 

-0.0900 

-0. 1 

3 




G„ 

-0.0871 

0.0203 

0.0667 

4 




q f (psi) 

0.633 

1.564 

1.745 

rad/sec) 

78.910 

74.646 

55.180 

GM (db) 

— 

0.798 

7. 993 

GM 3 (db) 

— 

0. 613 

8. 743 

GM (db) 

— 

9. 231 

4. 148 

aq f /5G 2 

1.94 

6. 97 

-1.62 

dq/dG 

2. 81 

9. 30 

1.27 

5 V 5G 4 

5. 00 

45. 04 

-10.02 

3q /3G 2 

-3. 82 

-153.4 

34. 06 

3q f / 3G 3 

-15.86 

-39. 3 

12. 09 

1 ~_—- 


Table 3: Three Base-Line Design Cases with Different Control Gains 


AW is the added weight, AW is the allowed weight penalty and and A f are cost 
weighting parameters. The sensitivity derivatives of J are 


3J 

dp. 


\ W i 


A w (AW-AW j ) 3q f A f 0 _c I f ) £ aGM j (GM q -GMj ) 


- A 


f dp. 


-I 




( 22 ) 


j=2 


where p is either a structural variable p or a gain G. and w. is the weight 

s. 11 


penalty associated with p =1. 

s. 

i 


A steepest descent optimization algorithm has been used to perform the four 
optimization cases presented in Tables 4 and 5. Case I of Table 4 has been 
performed to find the optimal gains, limited by G=0.1, without changing the 
structure. The results are close to those of case 3 of Table 3. but with better 
q^, and gain margins. However, does not satisfy the required margins. The 
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parameter 

Case I 
gain only 

Case II 

Case III 

G 3=° 

Case IV 

cL_. 

0. 1 

0. 1 

0. 1 

0. 9 

AW(lb) 

0. 

0. 070 

0. 277 

0. 287 

q f (psi) 

2.031 

2. 191 

2. 183 

2. 620 

o)^( rad/sec) 

51.397 

59.307 

83.610 

66.420 

GM 2 (db) 

6.57 

6. 56 

9. 57 

6. 39 

GM 3 (db) 

7. 48 

6. 69 

— 

5. 24 

GM (db) 

7. 07 

8. 37 

8. 92 

8. 84 

G 2 

-0.0974 

-0.0974 

-0. 1 

-0. 09 

G 3 

-0.0914 

-0.0901 

0. 0 

-0. 09 


0.0481 

0.0435 

-0.0315 

0.0230 


Table 4: Four optimization cases. 


i 

zone 

ply 

direction 

Case II 

(°) 

Case III 

Case IV 

i 

1 

0 

0. 199 

1 . 

0 . 

2 

1 

45 

0.020 

1 . 

0 . 

3 

2 

0 

0.317 

1 . 

0 . 

4 

2 

45 

■ 0.015 

1 . 

0 . 

5 

3 

0 

0. 002 

0 . 

0 . 

6 

3 

45 

0. 098 

2. 

0 . 

7 

4 

0 

0 . 

0 . 

0 . 

8 

4 

45 

0. 157 

2. 

0 . 

9 

5 

0 

0. 014 

0 . 

0 . 

10 

5 

45 

0. 035 

2. 

0 . 

11 

6 

0 

0 . 

0 . 

0 . 

12 

6 

45 

0. 075 

0 . 

0 . 

13 

7 

0 

0 . 

0 . 

0 . 

14 

7 

45 

0 . 

0 . 

0 . 

15 

8 

0 

0 . 

0 . 

0 . 

16 

8 

45 

0 . 

0 . 

0 . 

17 

9 

0 

0. 538 

6 . 

4. 

18 

9 

45 

0 . 

0 . 

0 . 

19 

10 

0 

1.346 

0. 092 

4. 

20 

10 

45 

0 . 

0 . 

0 . 

21 

11 

0 

0 . 

0 . 

4. 

22 

11 

45 

0 . 

0 . 

0 . 


Table 5: Structural changes (p ) in Cases II to IV 

I 


simultaneous structure and gain optimization of Case II yields sufficient 
margins at the penalty of adding 0.07 lb to the structure. It can be noticed 
from Case II of Table 5 that the 0 plies in the wing root zones 9 and 10 and in 
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o 

the wing tip zones 1 and 2 are the most effective. The 45 plies of zones 3 and 
4 have also some effect while the others are either not effective or constrained 
to the minimum gage. Case III is the same as II but with G^O. The required 

margins are also met but with an additional weight penalty of 0.207 lb. The 
flutter frequency and the structural changes are substantially different than 
those of Case II. The 0 and 45° plies of zones 1 and 2, and the 45 plies of 
zones 3 to 5 achieve their maximum limits while all the others remain at or 
close to their minimum limits. Case IV is the same as Case II but with gain 
limits of G=0. 9. The optimization process stopped in this case when all the 
independent design variables reached their limits. The differences between the 
resulting structures of Table 5 demonstrate the importance of simultaneous 
structural and control optimization of aeroservoelastic systems. 

The variations of q , AW, p and p along the 10 step optimization 

1 S 6 S 21 

process of Case III are plotted in Fig. 6. This case has been performed with a 



Fig. 6. Variations of flutter dynamic pressure, added weight and two structural 
variables along the optimization path, Case II. 

34 state aeroservoelastic model (22 structural, 6 aerodynamic and 6 control 
states). The computation cpu time was about 5 min on a MicroVax 3200. The 
dashed lines of Fig. 6 are for the same optimization process but with all the 25 
data-base modes yielding a 62 state aeroservoelastic model, which took about 32 
cpu min. Figure 6 demonstrates the efficiency of the optimization process and 
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the accuracy of the reduced—siz© model relative to the full size one. 

Two NASTRAN model updates were performed with the p g values of Cases III 

and IV of Table 5. The entire aeroservoelastic modeling process of Fig. 1 was 
repeated for each case, followed by root-locus stability analysis. All the new 
stability characteristics were within 3% of the associated ones in Table 4, 
which demonstrates the accuracy of the presented approach. 

Conclusions 

A systematic method for data-base construction, modeling and structural 
and/or control optimization of aeroservoelastic systems has been presented. The 
method is applicable to any objective function that can be analyzed by the 
normal modes approach and analytically differentiated with respect to the design 
variables. Standard, commercially available finite-element and unsteady 

aerodynamic computer codes can be used to construct the aeroelastic data-base. 
Application of Minimum-State rational approximations of the aerodynamic matrices 
yields a low number of aerodynamic states relative to the number of structural 
states. The weighting algorithm used to select modes for model size-reduction, 
and the dynamic residualization method have been shown to yield high-accuracy, 
reduced-order models. Sensitivity derivatives of flutter dynamic pressure and 
control margins have been extended to deal with residualized systems. 
Expressions for an efficient calculation of aeroelastic effectiveness parameters 
and their sensitivity derivatives have been presented and used to implement 
aircraft performance requirements. The ensemble of all these analytical and 
computational techniques yields an efficient, high-accuracy optimization process 
which can be used as a practical design tool of realistic modern aeronautical 
composite structures and control systems. 
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